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ABSTRACT 

In the classical theory of accretion disks, Lynden-Bell & Pringle demonstrated that angular momentum transport 
and mass diffusion are associated with the generation of energy through viscous dissipation. In the seminal model 
of Chiang & Goldreich, the only heating mechanism for the outer regions of protostellar disks is assumed to be 
stellar irradiation. Here, we construct a new set of self-consistent analytical disk models by taking into account 
both sources of thermal energy. We analyze the non-isothermal structure of the disk across the mid-plane for 
optically thick disks, and use the standard two-temperature model in the case of optically thin disks. We deduce 
a set of general formula for the relationship between the mass accretion rate and the surface density profile. With 
this prescription, the evolution of the disk can be studied for a set of arbitrary initial and boundary conditions if 
desired. When a Minimum Solar Nebula surface density profile is prescribed, our results recover those of Chiang 
& Goldreich in the optically thin regions, but extend their work for the opaque regions of the disk. 

For the purpose of illustration, we apply our theory in this paper to determine the structure of protostellar disks 
around T Tauri stars under a state of steady accretion and derive the corresponding radial variation of disk proper- 
ties such as surface density and temperature near the mid-plane. Our model predicts that the overall photospheric 
temperature of the opaque disk regions is typically smaller than previously thought. This modification is counter- 
intuitively related to a reduced flaring angle of the disk when viscous heating near the mid-plane is taken into 
account. We calculate the position of the snow line around a sun-like T Tauri star, and deduce that it can evolve 
from well outside 10 AU during FU Orionis outbursts, to ~ 2 AU during the quasi-steady accretion phase, to 
the present-day orbital radius of Venus when the accretion rate falls to about 10~ 9 M Q yr" 1 , and finally re-expand 
slightly beyond 2.2 AU during the protostellar- to-debris disk transition. This non-monotonous evolution of the 
snow line may provide some novel and deterministic explanation for the origin and isotopic composition of water 
on both Venus and the Earth. In addition, in classical T Tauri stars we infer the presence of a marginally opaque 
region between the optically thin and thick regions, with a surface density profile varying as r"V 2 , as in a Min- 
imum Solar Nebula model. The temperature is both vertically and radially constant and up to 40% higher than 
that in the region immediately within. This transition may lead to an upturn in the SEDs in the MIR (24-70/im) 
wavelength range. In the optically thin, outermost regions of the disk we find that the surface density profile of 
the dust varies roughly as r~ l , which is consistent with mm observations of spatially resolved disk of Mundy et al. 
(2000) 

Subject headings: accretion disks - methods: analytical - solar system: formation 



1. INTRODUCTION 

The coplanar orbits of the planets around the Sun inspired 
Laplace (1796) to postulate the nebula hypothesis for solar sys- 
tem formation. The first evidence for these long anticipated 
protostellar disks was inferred from continuum IR emission at 
a level well above the stellar black-body spectrum for a large 
fraction of stars in young clusters (Adams et al. 1987). Some 
flattened structures have also been resolved with mm imaging 
(Beckwith & Sargent 1992). The IR and mm excess emissions 
are attributed to the reprocessing of stellar light by dust grains. 
Since these disks are assumed to be the cradle of planet forma- 
tion, the determination of their structure and evolution has been 
a central task in the quest to reconstruct the origin of planetary 
systems. 

Obtaining accurate models of the thermodynamical structure 
of protostellar disks is important both from observational and 
theoretical points of view. Firstly, while spatially resolved ob- 
servations are possible for disks around nearby stars, the ma- 



jority of the available data is in the form of spectral energy 
distributions (SEDs), and can only be fully interpreted with a 
thorough understanding of the underlying disk properties. Sec- 
ondly, the details of the interaction between the gas and the con- 
densed heavy elements (in any form from dust grains to plan- 
ets) depend sensitively on the temperature and pressure gradi- 
ent of the gas within the nebula. In short, to understand both 
observationally and theoretically the processes associated with 
planetary formation requires first and foremost a better under- 
standing of the basic disk structure. This is the task we set out 
to do in this paper. 

1.1. Disk SEDs and grain properties 

In recent years, multi-wavelength observations have pro- 
vided a rich data bank of SEDs for most of the young stellar 
objects within several nearby clusters (cf Lada et al. 2006). 
Detailed numerical modeling in terms of reprocessed stellar ir- 
radiation has been successful in reproducing observations (Gol- 
dreich & Chiang 1997, Bell 1999, D'Alessio et al. 2001, Cal- 
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vet et al. 2002, Dullemond & Dominik 2004). These analyses 
show that SEDs are particularly sensitive to three properties of 
the disks: 1) the nature of the dust grains (density, size dis- 
tribution, physical and chemical composition), 2) the thermal 
and chemical structure of the gas, and 3) the global geometrical 
structure of the disk such as the inclination, degree of flaring, 
size, presence/absence of a hole (Dullemond et al. 2006). 

Current models of dust particles in disks usually assume that 
they have the same power-law size distribution as that of the in- 
terstellar grains (Mathis, Rumpl & Nordsieck 1977, MRN here- 
after), but with varying upper and lower cutoffs. This power law 
of the MRN distribution implies that the total cross sectional 
area is dominated by the contribution from the smallest grains 
(below micron size) whereas the large grains (> 1 mm) contain 
most of the solid-phase mass of heavy elements. Consequently 
the disk SEDs in the wavelength range observed by Spitzer are 
most sensitive to the small grains, which are tightly coupled to 
the gas both thermally and dynamically. 

The objective of this paper is to construct a set of compre- 
hensive analytical models of the thermodynamical structure of 
the disks, which can easily be utilized for the statistical analysis 
of large SED datasets and a more physical interpretation of the 
planet formation environment. 

1.2. Gas-dust interaction 

In parallel with the observational aspects of the study of 
the disk SEDs, large theoretical efforts have recently been de- 
ployed toward understanding the growth of small grains as well 
as the motion of the larger particles (and small planetesimals) 
within the primordial solar nebula (Weidenschilling 1977, Wei- 
denschilling & Cuzzi 1993). The properties of the dust grains 
have recently been thoroughly studied since evidence for dust 
growth (D'Alessio et al. 2001, Shuping et al. 2003, McCabe et 
al. 2003) or local (spatial) accumulation (Greaves et al. 2005) 
could be the first signs of ongoing planet formation. 

While the tiniest dust grains are dynamically very strongly 
coupled with the gas and largely follow its evolution, interme- 
diate size objects ("pebbles" to "boulders") have a different in- 
trinsic motion which is sensitively dependent on the structure 
(temperature, pressure gradient) and dynamical nature (turbu- 
lent versus laminar) of the nebular gas (Morfill & Voelk 1984, 
Supulver & Lin 2000, Youdin & Chiang, 2004, Ciesla & Cuzzi 
2006). 

It is thought that gas drag causes m-size objects to spiral into 
the central star within a few hundred orbital periods (Adachi 
et al. 1976), but the ubiquitous presence of planets around at 
least 10% of all observed stars (Marcy et al. 2000) suggests 
the existence of a reliable, yet still unidentified mechanism for 
the retention of at least a few percent of the total content of 
heavy-elements passing through the nebula. Meanwhile, up- 
per limits on this retention efficiency have independently been 
set through the observed homogeneity of the heavy-element 
abundances among sun-like stars in older clusters such as the 
Pleiades (Wilden et al. 2002), the Hyades (Quillen 2002) and 
IC4665 (Shen et al. 2005). These studies suggest that the to- 
tal mass of heavy elements retained by the protostellar disks 
towards planet formation is at most a few times that contained 
in a Minimum Solar Nebula (MSN, hereafter; Hayashi et al. 
1985). 

A self-consistent model of the thermal structure of proto- 
stellar disks is essential to understand these rather strongly 
constrained estimates of the heavy-element retention efficiency 



during planet formation. 

1 .3. Protostellar disk models 

The thermal structure of the gas is determined by a bal- 
ance between the generation of heat (through stellar irradiation 
(Adams et al. 1987) and viscous dissipation (Lynden-Bell & 
Pringle 1974), its transport (through radiative diffusion and ad- 
vective transport) and finally cooling (through emission by the 
gas itself, or more importantly by the dust grains with which the 
gas is thermally connected). In existing analytical models of 
protostellar disks, the energy source is usually either attributed 
to viscous dissipation (Lin & Papaloizou 1980, 1985, Ruden 
& Lin 1986) or stellar irradiation (Adams et al. 1987, Chiang 
& Goldreich 1997). More recently, numerical models of the 
disks have begun to take both sources into account (Bell 1999, 
d'Alessio et al. 2001, Dullemond & Dominik, 2004). 

In this paper, we propose a new simple kind of analytical 
model to calculate the two-dimensional (r,z) thermal structure 
of the disk taking into account both energy sources. The na- 
ture of the model enables us to derive analytical formula for all 
of the disk quantities, and derive scaling relations between the 
disk properties and the host star properties. 

Our work is largely inspired from the seminal paper by Chi- 
ang & Goldreich (1997) (CG97 thereafter). As they suggest, 
most of the stellar radiation is intercepted by a superheated 
layer of grains which re-emits part of it towards the disk mid 
plane. We also adopt an isothermal structure in the optically 
thin regions above the disk photosphere (in the case where the 
mid-plane is optically thick) and in the marginally opaque or 
optically thin outer regions of the disk. However, we introduce 
two major modifications to their work. 

Firstly, we model the optically thick regions in more detail 
by considering their non-isothermal structure. A warmer mid- 
plane is a natural consequence of heating by viscous dissipa- 
tion (in the case of optically thick disks) combined with down- 
gradient transport towards the photosphere. We assume a poly- 
tropic structure for the disk stratification in the direction across 
the mid-plane. This assumption can be shown to be consistent 
with some possible heat generation and transport mechanisms, 
as for instance in the case of viscously heated convectively un- 
stable disks. It can also naturally be considered an approxima- 
tion to the true stratification of the disk, with a polytropic index 
to be chosen from best fits of numerical simulations. 

Radiative heating by the central star is well-known to be the 
dominant heat source throughout most disks; nonetheless for 
completeness we also consider regions of the disk in which the 
dominant heating process is through viscous heating, which 
naturally occurs close to the central star. While this effect 
is not particularly important for classical T Tauri stars (M < 
10~ 7 M Q /yr) outward of 1AU, the viscously dominated region 
can extend up to tens of AU for much larger mass accretion 
rates. 

Secondly, while CG97 specifically assumed an MSN surface 

density profile (Smsn = 1000r^u 2 g/cm 2 ) we derive instead self- 
consistent formula between the disk structure and the mass ac- 
cretion rate. The evolution of the disk can then be studied for a 
set of arbitrary initial and boundary conditions if desired. In this 
paper we restrict our study to the case of steady accretion (with 
a radially constant M) for the sake of simplicity. While still 
being an idealized assumption, this is the state towards which 
accretion disks relax should a constant supply of material be 
provided near their outer rim. Consequently, it is a better ap- 
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proximation to the disk structure for the inner regions of the 
disk (and presumably out to a few hundred AU) (Lin & Pa- 
paloizou 1980, 1985) than the standard MSN power law surface 
density profile usually used. 

1.4. The importance of the snow line 

The temperature of a typical T Tauri disk decreases from over 
1000K within a fraction of 1AU to well below 100K beyond 
10AU (see §6). In the low pressure and tenuous environment 
of protostellar disks, ice grains sublimate between 150K and 
170K. The snow line is generally defined to be the location 
where the gas temperature attains these values (Hayashi et al. 
1985). 

As an example of the potential of the new disk model we 
propose, we present new results on the evolution of the snow 
line around solar-type T Tauri stars. This analysis is much 
more than an academic exercise, however, since the snow line 
is thought to play a very important role in the formation of plan- 
ets. 

The presence of a snow line somewhere in the disk sep- 
arates regions primarily containing refractory silicates from 
those containing volatile ices. Since the relative abundances 
of the lighter elements is at least a factor of four larger than 
that of metals, the supply of grains as planet building materi- 
als is much richer beyond the snow line. The large difference 
between the typical radial velocity of heavy-elements in con- 
densed form and in gaseous form suggests the possibility of 
accumulation of material near the snow line (Takeuchi & Lin 
2005, Ciesla & Cuzzi 2006). The thermal "anomaly" in the 
disk structure caused by the presence of the snow line may also 
be associated with a flattening of the pressure gradient, which 
could slow down the radial migration of the larger particles, and 
further enhance the local accumulation of material. Finally, this 
boundary may separate the domain of terrestrial planet versus 
gas giant planet formation (Ida & Lin 2004). 

Observational evidence for the presence of a snow line in 
disks around other stars is scarce. Measurements of the loca- 
tion of the snow line in the protosolar nebula have been inferred 
from the actual water content and suggestions of aqueous al- 
teration in various families of meteorites (Lunine 2005) to be 
between 2.5 and 3 AU, within the asteroid belt. 

We will show that while our model naturally reproduces its 
current location, it also predicts large excursions in the posi- 
tion and extent of the snow line during the early stages of the 
protosolar nebula. We will discuss the consequences of these 
findings for the formation and hydration of planets in our solar 
system. 

1.5. Outline of the paper 

The paper is organized as follows. In §2, we briefly recapit- 
ulate the total energy balance in the disk as a function of the 
dominant heating and cooling processes. We largely follow the 
footsteps of CG97 except where explicitly mentioned. In par- 
ticular, we use their prescription for the absorption and emis- 
sion properties of the dust grains, which are both sensitively 
dependent on the typical wavelength of the photons considered. 
This leads to a complex vertical and radial structure of the disk, 
which is also presented in §2. 

In §3 and §4, we derive analytical formula for the density 
and temperature profile in each of the layers for each of the re- 
gions considered, and deduce the two-dimensional structure of 
the disk. The methodology followed consists in calculating the 



position of the superheated dust grains as a function of the total 
optical depth of the disk, deducing the flaring angle of the disk 
and the associated intercepted stellar flux and equating it with 
the energy lost through radiation to deduce the thermal struc- 
ture of the disk. Once obtained, we can then evaluate the mass 
accretion rate. In the steady-state accretion approximation, the 
constant mass accretion rate assumed is then used to calculate 
the total surface density of the disk, and to deduce the total 
optical depth which closes the system of equations. Naturally, 
some of these steps are simplified in the isothermal case, and in 
the case where the dominant energy source is viscous heating. 
By following this method, we ensure that the disk structure is 
determined in a self-consistent manner. 

In §5, we estimate analytically the location of various crit- 
ical radii as a function of the accretion rate and the mass of 
the central star. These radii include those separating the radia- 
tively dominated optically thin, marginally opaque, optically 
thick and finally viscously dominated regions of the disk. We 
also obtain a rough estimate of the radial location of the snow 
line from the radius below which the disk mid-plane tempera- 
ture exceeds 160K. 

In §6, we apply our prescription to compute the structure of 
disks with a range of accretion rates around T Tauri stars. We 
confirm that throughout most of the disk, the flaring angle in a 
steady-state accretion state increases with radius. In the opti- 
cally thin region, we recover the scalings of CG97 for the tem- 
perature and scale height profiles, since these are independent 
of mass accretion rate. However, our surface density profile 
falls off with radius much less steeply (as r" 1 ) than that assumed 
for the MSN. This finding is consistent with recent observations 
of spatially resolved disks (Mundy et al. 2000, Looney et al. 
2003). In the marginally opaque region, we also coincidentally 
recover the scalings of CG97 since our model predicts in this 
region a surface density profile varying as that of the MSN. The 
temperature is interestingly found to be significantly larger than 
in the optically thick region just within: thus we predict the ex- 
istence of a "hot ring" at radii ^20-30 AU around classical T 
Tauri stars. In the optically thick regions our model predicts a 
significantly flatter disk, with notable consequences on the sur- 
face density, photospheric temperature profile and on the posi- 
tion of the snow line. We deduce three important results: the 
snow line can be very extended, it has a largely non-monotonic 
evolution, and can move (for a solar-type T Tauri star) as far in 
as the current orbit of Venus. 

Finally, we discuss the assumptions of the model in §7, as 
well as its implications for planetary formation and observa- 
tions of protostellar disks. 

2. BASIC MODEL ASSUMPTIONS AND METHODOLOGY 

2.1. Basic framework 

Our model follows a standard 1+1D methodology which im- 
plicitly assumes that the disk is thin, with independent vertical 
(i.e. normal to the disk plane) and radial structures. The dust, 
which largely governs the thermal structure, is initially assumed 
to be fully mixed with the gas; this assumption will be relaxed 
in future work. In the present context, we adopt a prescription 
in which 

Pd=Zp, (1) 

where p is the gas density, pd is the dust density, and the "metal- 
licity" Z is in principle a slowly varying function of radius 
and time. However, note that the magnitude of Z may change 
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abruptly near transition fronts (eg, it decreases by a factor of 
four as ice grains cross the sublimation line). We will consider 
the details of the disk structure near the snow line in a follow- 
up paper. In addition, note that pd refers to the heavy-element 
mass density contained in small particles only (with settling 
times which are much larger than the turbulent diffusion time); 
there could be a significant fraction of "undetectable" heavy- 
elements condensed in larger-size objects, that the model pre- 
sented here neglects. 

The disk has two principal thermal energy sources: radiative 
heating (from the central star) and viscous heating (from turbu- 
lent dissipation associated with momentum transport). In inter- 
mediate and outer regions of typical disks around T Tauri stars, 
radiative heating is by far the largest contribution to the total 
energy budget and essentially determines the thermal proper- 
ties of both the grains and the gas in the disk. In optically thick 
regions, even when radiative heating dominates the total energy 
budget, viscous heating plays an important role in maintaining 
a warmer mid-plane layer, and determines the vertical stratifica- 
tion of the disk. Finally, much closer to the central star, viscous 
dissipation becomes the main agent in heating the disk. 

The effectiveness of the conductive thermal coupling be- 
tween the gas and the dust is crucial to the determination of 
the disk structure. This coupling depends essentially on the 
wavelengths considered, as well as on the collision time-scale 
between the dust particles and the gas molecules. CG97 dis- 
cussed these processes in their seminal work; we will not repeat 
them here. 

For computational convenience and reasonable accuracy, we 
assume as they do that the emissivity and opacity of the grains 
interacting a black-body photon spectrum peaked at tempera- 
ture T; are 



and Ki = Ky 



(2) 



When applying the model to T Tauri stars (see §6), we shall 
adopt the value (3=1 although all the formula are given for ar- 
bitrary values of (3. In the same spirit, we shall also adopt the 
approximate value of k v = lcm 2 /g suggested by the work of 
d'Alessioef al. (2001). 



dust grains, the H2 molecules, and the coolant molecules (such 
as H2O, and others). However, the gas temperature in layer 1 
has relatively little impact on the vertical stratification of the 
lower regions of the disk, and for simplicity we shall set it to 
either the temperature near the mid-plane T m (in the optically 
thin case), or the effective temperature of the disk T e (in the 
optically thick case). 

Super-heated dust grains re-radiate half of the absorbed flux 
towards infinity and half towards the disk mid-plane. Where 
this radiation is absorbed depends on the optical depth of the 
mid-plane layer to re-processed IR photons. Thus we con- 
struct layers 2 and 3 where layer 2 is optically thin to the ra- 
diation emitted by the superheated dust-grains, while layer 3 
is optically thick to it (note that layer 3 only exists for large 
enough surface density). The two layers are separated by the 
dashed line, defined as the position where the optical depth of 
the grains to the radiation from the superheated layer reaches 
unity on a path perpendicular to the midplane. 

Finally, the heated mid-plane grains re-emit towards infinity 
at some temperature lower than T s still. Again, the total outgo- 
ing flux depends on the optical depth of the mid-plane to its own 
emission: if the surface density is large enough, the disk is op- 
tically thick to its own emission and we construct a fourth layer 
(layer 4) delimited by the photosphere (solid line). At the pho- 
tosphere, the disk emits a black-body spectrum with tempera- 
ture r e which is always lower than T s . Owing to the difference 
between T e and T s there exist a small difference between the 
mid-plane layer photosphere (solid line) and the edge of layers 
2 and 3 (dashed line). While this difference is usually negligible 
(and neglected) when the inner disk is optically thick, there can 
exist an intermediate regime (also discussed by CG97) where 
the mid-plane layer is optically thin to its own radiation but op- 
tically thick to the radiation from the superheated layer; in this 
case, layer 4 vanishes. We refer to this region as the marginally 
opaque region. When the mid-plane layer is also optically thin 
to the superheated grains, layer 3 vanishes and we refer to this 
region as the optically thin region. Note that for even smaller 
total dust content, the whole disk becomes optically thin to the 
radiation from the central star. We do not model this last phase 
in this paper. 



2.2. Layered structure of the disk 

Following the model of CG97, we identify several layers in 
a typical disk. The detailed structure of these layers in various 
disk regions is schematically illustrated in Figure 1 . 

The top layer of the disk is referred to as layer 1 where the 
grains are directly exposed to the stellar photons. Since sub- 
mm size dust grains are inefficient radiators at IR wavelengths, 
and dust-gas collisions are too scarce to account for significant 
cooling, grains in layer 1 are "super-heated" to a temperature T s 
satisfying roughly 
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(3) 



where L + is the stellar luminosity, r is the distance from the cen- 
tral star, T s is the temperature of the superheated dust grains and 
e s is their emissivity at that temperature. 

The energy balance in equation (3) is valid as long as the 
grains are directly heated by the central star. Thus, layer 1 ex- 
tends from the central star outward until a visual optical depth 
1 is reached (dotted line). The temperature of the gas in this 
layer is difficult to determine without a complete description of 
the (dynamical, thermal, and chemical) interaction between the 



2.3. Variation of the layered structure with radius 

The total energy balance in the various radial regions of the 
disk is now described in more detail. 

Optically-thin disk regions. 

In the low-surface density regions, mostly at large disk radii, 
the entire mid-plane layer is optically thin to the radiation from 
the superheated grains as well as to its own emission. In these 
regions, layers 3 and 4 vanish while layer 2 is preserved. The 
grains in the mid-plane (now layer 2) are heated by reprocessed 
radiation from the superheated grains located on both sides of 
the disk. As a result, both gas and individual dust grains have a 
temperature T m given by 



A s 4 
^e s • — - — ~ = e m cr I , 



(4) 



s 2 ~4irr 2 

where A s is the total emitting area-filling factor of superheated 
grains above z s (i.e. the fraction of solid angle covered by 
grains), T m is the mid-plane temperature (the gas is assumed to 
be isothermal in this region), and e m is the corresponding emis- 
sivity. Note that the factor of e s on the LHS accounts for the 
fact that the grains are inefficient absorbers at IR wavelengths. 
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FIG. 1 . — Cartoon representation of the layered model for the three regions considered. From left to right, the optically thick, marginally opaque and optically 
thin cases. In all regions, the top layer (1) contains superheated grains directly exposed to the (mostly visible) radiation from the central star. Its lower boundary 
(the dotted line) is referred to as z s (r) and is calculated by setting the optical depth from the central star to the line to unity (see the Appendix). The superheated 
grains re-radiate in the IR both towards infinity and towards the mid-plane. Some of the mid-plane grains are directly irradiated by the re-processed light from the 
superheated grains (layer 2). In the optically thick and marginal cases however, part of the mid-plane (layer 3) is optically thick to this emission. Grains in layers 2 
and/or 3 in turn re-process the absorbed light and re-emit it at slightly longer wavelengths still towards infinity. When the disk is optically thick to its own emitted 
IR light a fourth layer is constructed, and the emission occurs from the photosphere (solid line) located at Ze(r). Note that unless the fourth layer is entirely absent, 
layer 3 is very thin and can be neglected. 



We also include the additional factor of 2 accounting for the 
emission from both layers, which was neglected by CG97. 

Marginally opaque disk regions. 

In an intermediate range of surface density, the mid-plane layer 
is optically thick to the radiation from the superheated dust 
grains, but optically thin to its own emission. In this region, 
the disk is again assumed to be isothermal and layer 4 does not 
exist. Taking into account for the energy budget for the entire 
disk, the grain and gas temperature are determined by 



(5) 



where A m is the total area filling factor of grains in the mid- 
plane layer. 

Opaque disk regions. 

In the limit where the surface density is sufficiently large for the 
mid-plane layer to be optically thick to its own radiation, T m 
cannot be directly determined from the stellar flux. The total 
energy budget determines instead the effective temperature T e 
since the photosphere is reasonably close to the location where 
the reprocessed radiation emitted by the superheated grains is 
deposited. Thus 
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where T e is the temperature (of both dust and gas) at the pho- 
tosphere. Note that this equation assumes that the contribution 
from viscous heating is negligible, which can be shown to be 
valid outward of about 1 AU for a typical T Tauri star. If we are 
interested only in regions outward of the snow-line anyway, this 
approximation is justified. However, in the innermost region 
of the disk where viscous dissipation is the dominant heating 
mechanism the total energy budget becomes 



-Mnl = 2aT e \ 



(7) 



instead, where M is the local mass accretion rate, and f2 K is 
the local Keplerian angular velocity. This equation relates the 
total energy generated by viscous dissipation within the disk 
mid-plane to the emitted flux on both sides of the disk. 

2.4. Radial power law: notations 



For simplicity, in each of the four regions, we write the mid- 
plane quantities as 

T m (r) = Tr%j, (8) 
Pm(r) = pr 1 ^ . (9) 
where here and all that follows, we set r^v = r/(lAU). This 
uniquely defines the indices p and q, in a manner consistent 
with the notation of Takeuchi & Lin (2002). The resulting 
isothermal disk scale height is therefore obtained as 

_ 2+3 

h(r) = hr^. (10) 
Note that q > -1 corresponds to a flaring disk, whereas q < — 1 
corresponds to a self-shadowed disk. In all that follows, we 
assume that the disk is flaring and check this assumption for 
self-consistency. The disk surface density profile is 



£ = 2>' 



AU • 



(ID 

For instance, one can set s = 3/2 to recover the MSN model. For 
a steadily accreting disk on the other hand we determine S from 
the requirement that the total mass accretion rate M = 37Tf t £ 
be constant with radius, where we have adopted a standard a- 
prescription for the eddy-viscosity v t with value 

u t = a t c m h. (12) 

In that case 

s = q+l, (13) 



— _ M (1AU) 3/2 



(14) 



3. STRUCTURE OF OPTICALLY THIN AND MARGINAL DISK 
REGIONS 

In the outer regions of the disk, the total surface density of 
grains is sufficiently low for the disk to be optically thin to its 
own radiation. In this case, layer 4 vanishes. We assume, as in 
CG97, that the gas is isothermal in the direction normal to the 
plane of the disk and has a temperature T m . 

Integrating the equation for hydrostatic equilibrium yields 
the standard Gaussian vertical density distribution for the gas: 

8 , (15) 



P(r,z) = p m (r)cxp 



with h 2 = 



2h 2 (r)J "' 7 fi2 M 0| 
where c m is the sound speed for a gas at temperature T m , with 
mean molecular weight ^ and adiabatic index 7; R g is the gas 
constant. When the dust and the gas are fully mixed, the in- 
terface z s between layer 1 and layer 2 (which can be thought 
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of as the position of the superheated dust layer) is obtained by 
solving the equation 



f Zp(r',z')Ky6l = l 
Jo 



(16) 



on a straight path dl from the central star to the point (r,z s ) 
(namely, z' = (z s / r)r'). Here, Z is the assumed metallicity Z nor- 
malized to the value used for the determination of Ky, namely 
1%. Thus Z = Z/0.01. Note that it can also be regarded as 
~ 10 [Fe / H1 where [Fe/H] is the metallicity normalized with re- 
spect to that of the Sun. Assuming the disk is flaring, one 
can approximate this equation using the variant of a stationary 
phase method (see Appendix), which yields the equation (valid 
at every radius r within the region) 



Zp m liKy exp 



2h 2 



h L r dlnr 



where we have defined a as the grazing angle between the in- 
coming radiation from the star and the local stratification: 



d_/h 
dr \ r 



(18) 



Note that this approximation is only strictly valid in the limit 
where z s is larger than a few scale heights. By inspection, one 
can note that the second term on the RHS of equation (17) 
is negligible within the same limit unless the disk has a sud- 
den jump in gas density or metallicity (this could occur near a 
planet-induced gap, or near the successive sublimation zones). 
Deferring the study of these transition zones to a later paper, we 
find that the ratio z s /h is otherwise given by 



Zs 

h 



2LambertW ^ 



2ira 



1/2 



(19) 



where Ty ot = ZEkv/2 is the total visual optical depth from the 
mid-plane to infinity on a path perpendicular to the disk. Given 
that the disk is isothermal, E = \/2~Trp m h. 

Note that the LambertW function is defined as the inverse of 
the function f(x) = xexp(x), or in other words, if y = xexp(x) 
then x = Lambert W(y). For ease of computation, and for large 
values of y (for y between 10 and 10 5 say), a fairly good ap- 
proximation (to within about 10%) of the LambertW function 
is given by 

LambertW(y)~0.81n(y) . (20) 

This can be seen in Figure 2. 

We therefore note that as long as Ty 0t 
is a slowly varying function of radius with 



2na the ratio z s /h 



= C(r): 



3.7 log 



10 



lira 



1/2 



(21) 



Calculating C{r) exactly in the case of a classical T Tauri star 
reveals that C varies between 3 and 4 in the outer regions of the 
disk, confirming the ansatz (z s ~ Ah) of CG97. However, for 
other stellar types or for other mass accretion rates, we find that 
the superheated dust layer position varies widely (with C taking 
any value between 1 and 6). 

Once the magnitude of z s is determined, we need to identify 
the surface area filling factor A s of the reprocessing material 
contained within the superheated layer which is equal to the 
visual optical depth r v . This is done through the integral 



f 



ZpKydz = T v 



(22) 



Using the same stationary phase approximation (see Ap- 
pendix), we evaluate the integral to get 

( z 2 \ h 

A s ~ Zp m Kyhex V ( -^2 J - • (23) 
and combining this with equation (17) yields 

(24) 

h 

Note that this expression differs from the one suggested by 
CG97 by a geometrical factor of z s /h = C(r); this difference 
arises in choosing to determine A s from a path that is perpen- 
dicular to the mid-plane rather than perpendicular to the disk 
stratification. We prefer the former interpretation as more con- 
sistent with the 1+1D assumption, although the latter is also 
acceptable within the degree of approximation associated with 
this analytical approach. 

3.1. Optically thin disk regions 

To compute the radial structure of the disk, we combine 
equations (4), (24), and (18), and ignore the slow dependence 
in r of C(r). This derivation yields an equation for h(r)/r for 
instance, of the general form described and solved in the ap- 
pendix: 

which can be integrated to yield equation (10) with 



h_ 
AU 



1 2/3 2 + 15/3+28 / p GM 



C /3 2 + 4/3+8 \R g R+Ti, 



4+li 



7+2/3 



V1AU 



__8+4/3+/3j_ 
28+15,3+2,3 2 



q + 3 P 2 +4P+S 



2/3 2 + 15/3+28 



+ 1 



(26) 



If /3 = 1 then q = -19/45, with h oc r^ 45 and the flaring an- 
gle a oc r 1 ^ 45 . This result is independent of the surface density 
distribution of the gas, and is similar (within factors of order 
unity) to the expression derived by CG97. If we require in ad- 
dition that the disk be in a state of steady accretion, the surface 
density distribution can be deduced directly from (11), and in 
that case s = 97/90 ~ 1. This result is particularly interesting 
in the light of the observational results of Mundy et al. (2000), 
who found that the surface density profile of spatially resolved 
disks varied as r~^ v rather than r~^ 2 . 

This analysis is only valid in the region delimited from the 
inside by the inequality r' ot < 2/3 (where T^ ot = ZEk s /2), and 
from the outside by the inequality Ty 1 3> a (or h ~ r, whichever 
happens closest to the central star). Note that as the total optical 
depth drops below a, the superheated layer becomes progres- 
sively more and more optically thin to the radiation from the 
central star. This is expressed mathematically by C — > 0. When 
this happens, all the grains directly exposed become "super- 
heated". The temperature of the gas, however, is more difficult 
to evaluate because of the long collision time scale between 
grains and gas molecules; as a result, it is unclear whether 
the disk retains its flaring structure or becomes self-shadowed, 
although we suspect the latter is more likely. This transition 
would appear as an outer edge to the observable disk. 
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FIG. 2. — LambertW(x) function (solid line) and its approximation 0.81n(;t) (dashed line) for large values of x. Although this is not an asymptotic approximation, 
and is not uniformly valid for x — > oo, it is quite good in the range considered. 



In addition, when considering the case of steady-state accre- 
tion, it is important to bear in mind the underlying assumptions 
that (1) the disk has time to relax to the quasi-steady state, and 
(2) that a constant mass supply is indeed provided to the disk 
from the outer rim. Both these assumptions are difficult to jus- 
tify beyond a few hundred AU. 

3.2. Marginally opaque disk regions 

To compute the radial structure in the marginally opaque 
case, we need to determine the emitting surface of grains in the 
mid-plane region z £ [0,z s ] which is very well approximated 
by Ty ' in the limit where z s is larger than a few scale heights 
(which we proved to be self-consistently true): 

A m ^r = ^. (27) 

Using the same method as in the previous section (which in- 
volves using equations (5), (24) and (27)), we find that the disk 
scale height is 



h 



1AU 
q + 3 



2 7 + 2/3 ZSk v / li GM„ 
C2 + S + /3 2 \R~ g R^ 

9 + s + 3(3 



4+/3' 



7+2/? 



It.) 

iAuy 



2+13 
1-1 



(28) 



2 7 + 2/3 ' 
Interestingly, we find that the MSN scaling law is the one 

3 12 

consistent with steady accretion, with h oc r^ v and 

3 h 
s=- ,q = 0,a=- , (29) 
2 2r 

for any value of (3. This value of q implies that the disk is radi- 
ally isothermal in this region. In addition, S and h are obtained 
by solving simultaneously equations (28) and (14), so that 



1 AU 



MZk, x 



H GM* 

C 3ira ty /'yGMJt i , \R~ g R7^ 



4+/3' 



5+2/? 



— v 

1AU / 



and the value of the radially constant temperature is 

H GM* -2 

l m = ttH . 

Ro (1AU) 3 



(31) 



Again, this result recovers the expression obtained by CG97. 
The marginally opaque region terminates at the radius where 
r; ot = 2/3, where = Z Km S/2. 



4. VERTICAL STRUCTURE OF OPAQUE DISK REGIONS 

Disks are heated both by viscous dissipation and by stellar 
irradiation. When the mid-plane is opaque to its own thermal 
radiation, heat generated near the mid-plane can only be trans- 
ported towards the photosphere by turbulent or radiative down- 
gradient diffusion. This implies the emergence of a temperature 
gradient across the disk mid-plane, which must be calculated 
self-consistently. 

4. 1 . Thermal stratification in the optically thick layer 

Existing detailed disk models are constructed based on the 
numerical integration of this thermal equilibrium equation (cf 
d'Alessio et al. 2001; Dullemond & Dominik, 2004) 



'''' ( ''7j; (!! K> 



: V • (Frad + Fturb) . 



(32) 



Unfortunately the nonlinearities in the radiative contribution to 
the heat flux prevent simple analytical solutions. Instead, we 
adopt a polytropic structure for the disk which can either be 
deduced exactly assuming idealized heat transport mechanisms 
(such as effective convective transport leading to an adiabatic 
stratification), or considered an approximation to the solution 
of the full transport equation. 

Let's consider the standard polytropic equation 



(30) 



p = Kp 



(33) 
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where p is the gas pressure, n is the polytropic index and K is 
the polytropic constant. The equation of hydrostatic equilib- 
rium across the mid-plane is 

dp 



dz 



(34) 



where we have neglected the self-gravity of the disk. Combin 
ing these with a perfect gas equation of state 

_ R iP T 



/i 



(35) 



where R g is the gas constant and p, is the mean molecular 
weight, we obtain 

T{r,z) = T m {r)(\- 



p(r 1 z) = p m (r) 1 



H(rf 

z 2 
' H(r) 2 



(36) 



where H is uniquely related to the conventionally defined disk 
isothermal scale height 

A 

(37) 



(where 7 is the adiabatic index for the gas) as 
H 



— =6 = V2n + 2, 
h 



(38) 



Note that if the gas is composed of mostly molecular hydrogen 
with a solar composition, 

7= 1.4, fi = 2.4. (39) 

If we consider a convective polytrope for instance, we choose n 
such that 7 = 1 + 1/n or in other words, n = 2.5. 

Note also that H determines the location where, in principle, 
both the gas density and temperature would drop to zero. How- 
ever, this mid-plane solution is only valid within the optically 
thick interior of the disk, i.e within the interval z € [0,z e (r)], 
where z e (r) is the position of the photosphere. For the conve- 
nience of following discussions, we define the ratio 

A e = Zc/H . (40) 
Above the photosphere, the temperature and density profiles 
are matched onto an isothermal atmosphere solution with the 
photospheric temperature T e . The matching procedure involves 
solving the total energy equation to determine T e , which then 
uniquely yield the values of the mid-plane temperature and den- 
sity r m and p m . In order to do this, we must first determine the 
position of the photosphere, which can be done once the verti- 
cal stratification in the atmosphere is calculated. 

4.2. Disk atmosphere 

Above the photosphere, we assume the gas to be isothermal 
with temperature T e , and constant sound speed matched to the 
interface value: 

c 2 (r,z) = cl{r) = c 2 (r,z e ) = iM^l , (41) 

for z > z e - Integrating the equation for hydrostatic equilibrium 
we obtain the gas density: 

z 2 -A 2 H 2 " 



p = Pm (l-A 2 e ) exp - 



2hl 



where 



(l-A 2 ) = /z 2 (l-A 2 ). 



(42) 



(43) 



4.3. Location of the photosphere. 
To obtain z e , we must solve the equation 

2 

3 ' 



f 



Zpn e dz - 



(44) 



where K e = KyiTz/T*) 13 . Assuming that the dust is fully mixed 
with the gas (see equation (1)) yields the required equation for 

A e : 



(1"A>) 



2nk+,3+1/2 



erfc 



6A e 



7T Zp m hKy \ 



V^Af, 

-0 



exp 



2/)2 



Aid 

2-2A? 



(45) 



This expression for A e in terms of p m and T m at a given radius 
r seems rather complicated and requires in principle a numer- 
ical solution. However, there exist two limits in which a sim- 
ple analytical solution exists. In the strongly opaque (when the 
right-hand-side of equation (45) is very small) we expect z e to 
be located high above the mid-plane and very close to H (so 
that A e ~ 1). The nature of the polytropic/isothermal solution 
we are using implies that in this case, we expect 



1 - A 2 = e e < 1 so that A e ~ 1 ■ 



' 2 



(46) 



Substituting this ansatz into equation (45), and using the ap- 
proximation (76) we obtain 



/ (2/3)fl 2 /(n) 







(47) 



since in the limit where z e — H the total surface density of the 
gas can be approximated by 

/OO r-H 
p(r,z)dz~ / p(r,z)dz = 2p m (r)H(r)I(n), (48) 
00 J-H 



where 



/(«)= / (i-x 2 r 

Jo 



1 ( 1 
dx=-B -,„ + l 



(49) 



(B is a Beta function, see Gradshteyn & Rhyshik p. 343). From 
equation (47), we note that the approximation Ty 1 ^> 1 is con- 
sistent with £ e Cl. 

In the alternative limit, the right-hand-side of equation (45) 
approaches 1 (the disk is close to being optically thin), and we 
expect A e to be very small, in which case an approximate solu- 
tion can be obtained by linearizing (45) 



7T 1 

28 



1 



2/3 



-0 



(50) 



and using S ~ \J2-k p m h. 

The regimes of validity of the two approximations are illus- 
trated in Figure 3. For the purpose of computations, we will 
denote from now on the "strongly opaque regime" as the re- 
gion where we adopt the asymptotic approximation (47), and 
the "weakly opaque regime" as the region where we adopt (50) 
instead. Since we are not performing a rigorous matching be- 
tween the two regions, we artificially patch them instead at the 
radius where the RHS of (45) reaches 0. 1 . 
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FIG. 3. — Comparison between the true solution to equation (45) and the two asymptotic approximations. The x-axis represents the right-hand-side of equation 
(45). The solid line is the numerical solution, the dotted line is the approximation in the weakly opaque limit, while the dashed line is the approximation in the 
strongly opaque limit. 



4.4. Characteristics of the superheated dust layer 

The calculation of the position of the superheated layer is 
similar to the isothermal case, but with a different vertical den- 
sity distribution for the gas. As before, we use the stationary 
phase approximation assuming that the disk is flaring. In the 
weakly opaque limit, we retain the expressions derived in the 
optically thin case (namely equations (21) and (24)), whereas 
in the strongly opaque limit (see Appendix), we get 

z s ~H = 8h, (51) 

and as before 

A s = ^a ~ 9a . (52) 



^-a ~ 9a 
h 



It is interesting to note at this point that in a completely isother- 
mal model (as the one proposed by CG97), C(r) increases 
monotonically with Ty ot (equation (21)). Here instead we show 
that C(r) enters a different regime in the optically thick region, 
and has a constant value set uniquely by the gas properties con- 
tained in 9 (see equation (38)). 

4.5. Radial structure 

Weakly opaque region. Now that we have all the required 
information, we need to solve the total energy conservation 
equation. When the disk is nearly optically thin to its own 
radiation, the mid-plane temperature is very close to the pho- 
tospheric temperature: the disk is nearly isothermal (typically, 
e e ~ 1). In this case, we approximate, for simplicity, the total 
energy equation by 



2 r dr \r) r 2 \Tj ~ ^ 
Integrating this equation for the disk scale height yields 



1AU 



1/7 



"4/7 



(— 

V1AU 



-2/7 



q + 3 _ 9 
~~ 7 ' 



(53) 



(54) 



(55) 



so that q = -3/7 and for a steadily accreting disk, s = 15/14. 
This result, once more, naturally recovers the radial power law 
scaling found by CG97 since their work assumes that even op- 
tically thick disks are isothermal. Interestingly, we find that 
a = 2/7(h/r) which is notably smaller than the flaring angle in 
the nearby marginally opaque region (see equation (29). This 
difference indicates that the inclusion of the viscous dissipa- 
tion reduces the flaring angle compared to a purely radiatively 
heated disk. This result is not surprising. The scaleheight of a 
standard disk is determined by the balance between the pinch- 
ing force of gravity, and the vertical pressure gradient which is 
related to the local temperature. Thus the flaring angle is re- 
lated to the radial variation of the ratio of these two quantities. 
In the case of a radially isothermal disk, the flaring angle is 
maximal since gravity quickly decreases outward. In the case 
of radiatively heated disks, the addition of a negative temper- 
ature gradient implies that temperature increases inward, and 
reduces the flaring. Finally, the additional contribution from 
viscous heating to the mid-plane temperature further reduces 
the flaring. 

What is more counter-intuitive is that the inclusion of viscous 
dissipation can, in some regions, reduce the photospheric tem- 
perature of a disk compared to a model in which it is not taken 
into account: additional heating does not necessarily mean ap- 
parently hotter disks. This happens in regions where viscous 
heating is only important in maintaining a warmer mid-plane 
but does not contribute much to the overall energy budget. In 
this case, the reduced flaring angle implies a reduced inter- 
cepted stellar flux, and therefore a cooler photospheric temper- 
ature. 

Strongly opaque regions In what follows, we show how 
our formalism is able to extend the model of CG97 to non- 
isothermal, strongly opaque disks. In this new limit, we use 
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the fact that T e = e e T m with equation (47). After a few manipu- 
lations we get 



h 



1AU 



q+3 



2n + 2-As-209 
ln + l-[3 2 

4 6» 2 /(n)\ _7;; ^ 

3 ZSkv / 
9n + 9- 3/3-4* 



1AU 



fi GM* 

2«+2-2/3 



4n+4 
7I1+V-,) 



(56) 



(57) 



2 ln + 1-0 

In the case of a steadily accreting disk, we combine these 
equations with equation (14) to obtain 



h 



1 AU 



2n-2f3 9 
7n+15-/32 



s = 



15n+ 15-9/3 



2rc-2/3 
7n-£+15 



(58) 



(59) 



2(7n+15-/3) 

With /3 = 1, and n = 2, we obtain s ~ 0.65 and </ ~ -0.85, which 
confirms again that the disk is flaring. However, we note that in 
this case the flaring is even weaker (a ~ 0.1/z/r), and the disk 
surface is very nearly flat. Again, this geometry has very impor- 
tant consequences for the overall inner disk temperature, since 
a much lower fraction of the stellar luminosity is intercepted 
than in a model in which the disk is assumed to be isothermal 
(see the previous case). 

Viscously heated disks Finally, if the dominant contribution 
to the heating of the disk arises from viscous dissipation instead 
of stellar irradiation, manipulations of the total energy equation 
(7) imply that 



^3(< W 1/ 



1/4 



-3/4 



(60) 



which can then be used to deduce h, as usual: 

_ n+0+1 1_ 

8n+8 Ue 2 I{n)\ 2 " +2 / [iGM* 
V3kvSzJ \R g R*U 

, 3/3-n-l 

R, 



1AU 



3 Lace 

2~U 



-1/2 



1AU / 

q + 3 9n-30+9-4s 



(61) 



2 8(n+l) 

where we define here L acc = GM+M/R+. In the case of steady 
accretion, we can reduce this further to yield 

h ( 3 L acc \ 8 " +16 / 49 2 I(n) na t 



1AU 



acc 

2~u 



fx GM* 



V k v Z M 



T-y/jGM^kZ 



^R g R*T+ 
q + 3 9n-3/3+15 



1AU 



3 .'3-/?- 1 
8«+16 



(62) 



2 8n+16 
With the adopted fiducial values for n and /3, we find that the 
viscously heated regime is the only one that is not flaring: in- 
deed, in this case q = -1.125 and s = 0.375. However, flaring 
is not self-consistently required towards the determination of 
the disk structure in this region, so the formula given above 
is in fact valid. One may wonder, however, whether the vis- 
cously heated regions could shadow the innermost radiatively 



heated regions, as proposed by Dullemond & Dominik (2004) 
for HAeBe stars. We do not believe that this is the case for 
classical T Tauri stars: it should first be noted that the disk is 
still very nearly flat (q is very close to -1), and secondly that the 
radial extent of the viscously heated region is small (which is 
not the case for HAeBe stars). 



5. POSITION OF THE TRANSITION RADII 

In this section, we derive analytical expressions for the po- 
sition for the transition radii across the various regions for the 
specific case where the disk undergoes steady-state accretion. 
For simplicity we give the formula for the specific case where 
(3=1 and under the assumption of steady-state accretion, al- 
though they can also be easily derived in the general case. In 
practise, the transition radii are most easily computed numeri- 
cally. 



5.1. Transition from optically thin to marginally opaque 
This transition occurs when r t 5 0t = 1, or equivalently when 



Mr~ 



3Tra t y/^GMjl 2 2 \ r2 



= 1 , 



(63) 



where s and h are derived from equation (26). Solving for the 
transition radius r'^ m yields 



'AU 



/ ■ 



M 



2 3na lv ^GMJ^ 



1 45 / fj, GM+ 



C 13 \R g R+T* 



1AU 
(64) 

We see that, as expected, this transition moves inwards as the 
disks evolve and their mass accretion rate decreases. 



5.2. Transition from marginally opaque to optically thick 

This transition occurs when r™ t = Z~En m /2 reaches 2/3, or 
equivalently, when the RHS of (45) reaches unity. In a steadily 
accreting disk, with 0= 1, the right-hand-side of equation (45) 
reduces to 



2 2 1 (% 

3 V 7r Zp m hnv \ T-, 




(65) 



where K varies slowly between \J2-e (in the nearly-optically 
thin case) and 2I(n)9 (in the very optically thick case). This im- 
plies that for a given accretion parameters M and a t and given 
stellar parameters and M*, the transition radius at which the 
disk changes from optically thick to optically thin to its own 
radiation is entirely determined, and equal to 



'AU - 



Zk v M 



2 27ra t V ~/Rl R g T* 



2/3 



1AU 



(66) 



Again, as expected, this transition moves inwards as the disks 
evolve and their mass accretion rate decreases. It is interesting 
to note, however, that to a good degree of approximation both 
r t ^ m and r m ^ D are proportional to M 2 / 3 , which implies that the 
ratio r ^ m / r m ^ i s independent of mass accretion rate. 
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5.3. Position of the snow line 

The outer edge of the snow line is located where T m first 
reaches 160K. This transition can either happen in the optically 
thin region or in the optically thick region. In both cases, its 
position is given by the equation 



'AU 



/ 160# g (l AU) 3 



l/q 



(67) 



where h and q are region-dependent. In the optically thick case, 
(for (3=1, and n and 9 given by equation (38)) we find that 



'AU 



(68) 



so the snow line moves inwards as the mass accretion rate de- 
creases. Since this migration rate is slower than the rate at 
which the optically thin region expands, there must be a point 
where the snow line becomes optically thin. When this hap- 
pens, its location moves outward again to a radius which is a 
function only of the stellar parameters (with a weak dependence 
on the disk structure through C). When (3=1, 



'AU 



/ 160R g 



(l3c) 



ft GMA » R* 
R g R*Tj 1AU 



(69) 



The smallest possible radius attained by the snow line in the 
course of the disk evolution can be estimated simply by finding 
the position of the snowline while it is located in the weakly 
opaque region. We find that 



mm r 



AU 



: 0.6AU I — 



2/3 



-1/3 



— (70) 

Mr. ' 



This can be considered as the minimum radius at which a sig- 
nificant amount of vapor-phase water can be found in proto- 
planetary disks in the course of their evolution. Most impor- 
tantly, note that this value is independent of the value of k v se- 
lected, and of the nature of the polytropic solution selected: it is 
uniquely determined by the properties of the central star. How- 
ever, the value of M for which this minimal position is achieved 
does depend on these disk properties. 

In Section 7, we discuss in more detail the evolution of the 
position of the snow line as the disk evolves, and its conse- 
quences for planet formation. 

5.4. Transition from radiatively heated to viscously heated 

This transition occurs when viscous heating and radiative 
heating are of the same order of magnitude. In a steady-state 
accretion approximation, we deduce the transition radius from 
equating 

As L* 



—Mill 
4n 



2 4nr 2 



to be 



'AU 



9+3 



(71) 



(72) 



where the constant S is either 9 if the transition radius occurs in 
strongly opaque regions, or C if it is located in optically thin, 
marginally or weakly opaque regions. The exact values of h 
and q also depend on the radial region considered. 



6. EXAMPLES OF STAR/DISK SYSTEMS 

6.1. Examples of protostellar disks around a weak-line T Tauri 

star 

We now apply our model to dusty disks around T Tauri 
stars. We defer to later publications the discussion disks around 
the more massive Herbig AeBe stars and around the low-mass 
brown dwarves. Direct comparison between our model and the 
predictions of the CG97 model are presented in the discussion 
section. 

For the purpose of computing the structure of the disk, we 
select the following values for the global stellar and disk prop- 
erties (except when explicitly mentioned): 

U = \Lq,M* = \Mq,R*=R q , 
Z = 2 equivalent to Z = 0.02 , 
a, = 10 -3 , 

M = 2.4, 7 = 1.4, (73) 
n = 2, 

Ky = lcm 2 /g ,(3 = 1 

C = 4 , (74) 

and vary the mass accretion rate from M = 10~ 8 M Q /yr to M = 
10~ 10 M Q /yr, to cover a range of gas accretion rates from ob- 
served values for classical T Tauri stars (Hartmann 2001) down 
to low values presumably corresponding to a protostellar-debris 
disk transition. The results are shown in figure 4. 

Comparisons of the solutions for a particular value of M for 
two values of the polytropic index n and the grain opacity ny 
are shown in Figure 5. We note that while the choice of the 
polytropic index only has a small effect on the radial structure 
of the disk, the choice of the opacity naturally influences the 
transition between the optically thick and optically thin regions 
of the disk. 

For each accretion rate, three figures are shown: 1) the vari- 
ous relevant surfaces (h(r), z s (r) and z e (r)), 2) the relevant tem- 
perature profiles (T m (r), T s (r) and T e (r)) and 3) the surface den- 
sity profile. In each case, the superheated grains position and 
temperature are shown as a dotted line. Note that T s depends 
only on the nature of the central star: the profiles obtained in 
all three values of M for the particular T Tauri star chosen are 
identical. The position and temperature of the photosphere are 
shown as a solid line. The mid-plane temperature and its asso- 
ciated disk scale height h are shown as long-dashed lines. 

The three regions (from right to left, optically thin, 
marginally optically thin and optically thick - the latter sub- 
divided into weakly opaque, strongly opaque and viscously 
heated regions) are clearly visible in each figure. The transi- 
tion radii from optically thin to marginally opaque, and from 
marginal to optically thick are genuine physical transitions, and 
are marked by vertical dashed lines. We observe that the tran- 
sition from an optically thin to a marginal region is associated 
with a fairly continuous temperature profile, while that from 
marginal to optically thick has an important temperature dis- 
continuity. Note that in a real disk the temperature "jump" 
would be more diffuse (as a result of radial thermal and viscous 
diffusion processes unaccounted for in our 1+1D model), and 
could span a radial extent of up to one scale height. The arti- 
ficial transitions and associated discontinuities observed in the 
opaque region, (from left to right) between radiatively heated 
and viscously heated regions, and between the weakly opaque 
and strongly opaque regions, are mathematical artifact of our 
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FIG. 4. — Model prediction for a T Tauri star disk with accretion rates (from top to bottom) of M = 10~ 8 ,M = 10~ 9 , and M = 10~ 10 MQ/yr. The left-side figure 
shows the various important surfaces and heights, the central figure shows the relevant temperature profiles, and the right-side figure shows the surface density 
profiles. On the left, the dotted line marks the position of the superheated layer Zs(r), the long-dashed line marks the disk temperature scale height h(r), and the 
solid line marks the position of the photosphere z e (r). On the center, the dotted line marks the temperature of the superheated grains 7" s , the long-dashed line marks 
the temperature of the mid-plane T m and the solid line marks the effective temperature T e . On the right, the solid line marks our model prediction for E(r) while 

the thick dashed line is the standard MSN profile Emsn = 1000r A u 2 g/cm 2 . In all diagrams, the snow line for the mid-plane temperature appears as a vertical solid 
line. The vertical dashed lines (when present, from left to right) denote the physical transitions from the optically thick to the marginally optically opaque region, 
and then to the optically thin region. The vertical dotted lines marks the position of the artificial transitions between asymptotic regimes: near the central star (and 
in these diagrams, on the left of the snow line when present) the transition from viscously heated to radiatively heated, and near the marginally opaque region, the 
transition from the largely opaque region to the weakly opaque regions. 
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FIG. 5. — Comparison of the model predictions for the temperature profiles of a disk with M = lCT 8 MQ/yr around a solar-type star, for two values of Ky (from top 
to bottom Ky = 1cm 2 /g and k\ = 10cm 2 /g) and two values of the polytropic index n (from left to right n = 2 and n = 3). The four figures are otherwise generated 
using the same conventions as in Figure 4. Note that the choice of the polytropic index n has only a small qualitative influence on the solution in the optically thick 
regions of the disk (and naturally no influence on the optically thin regions), while the opacity k\ influences strongly the position of the transition between optically 
thin and optically thick regions of the disk. 
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crude patching procedure between various asymptotic approx- 
imations (in the former case, for the dominant heating mecha- 
nism and in the latter for equation (45)); they are represented 
by dotted vertical lines. One can estimate the true profile by 
smoothly matching one region to the next across these dotted 
lines by eye, or, should one desire a more accurate solution, 
actually solve the problem numerically. 

Finally, the vertical solid line represents the snow line for 
the mid-plane temperature. Across this line, the equations pre- 
sented here are only qualitatively valid, since we neglect all 
effects associated with the sublimation of the ices. We shall 
describe in detail a new snow-line model in a following pa- 
per. Note that the temperature of the superheated grains reaches 
160K at a far greater radius than the marked (mid-plane) snow 
line. The sublimation of the icy particles in the superheated 
layer results in a drop in Z above z s , but as shown in equation 
(21), this would only have a small effect on the value of C(r). 
Hence in what follows we ignore the effect of the successive 
sublimation lines on the structure of the superheated layer. 

Note that since the outer edge of the disk depends on a 
number of potentially important factors (initial radius, self- 
shadowing in the very optically thin case, total mass, thin-disk 
approximation, etc.), we have not marked it on the figures pre- 
sented. The reader should also bear in mind that the model has 
less predictive power in regions much beyond 100 AU, where 
the assumptions inherent to steady-state accretion are no longer 
satisfied; in these regions, the problem is best solved as an ini- 
tial value problem. 

Several notable qualitative features can be readily seen from 
these T Tauri star models and are also found in disks with a 
wide range of accretion rates, around stars with a wide range of 
stellar types. 

Firstly, we observe the presence of a wide, radially isother- 
mal region (corresponding to the marginally opaque region). 
Although the isothermal nature of this section of the disk was 
already inferred by CG97, we find that it also has the remark- 
able property of being notably hotter than the region immedi- 
ately within, so that the overall photospheric temperature pro- 
file is not a monotonically decreasing function of radius. The 
hotter temperature is a result of the larger flaring angle of this 
region compared to the one within, implying that a larger frac- 
tion of the stellar flux is absorbed locally. The presence of a 
"hot ring" could have very interesting consequences, both from 
an observational point of view and from the point of view of the 
dynamics of dust grains in the disk. This is discussed in detail 
in §7. 

Secondly, although a more rigorous snow-line model will be 
needed for more accurate quantitative discussions, it is already 
apparent in the models presented that the radial extent of the 
snow line can be considerable when it lies within the optically 
thick region. By radial extent, we imply the region between the 
radius at which the mid-plane temperature becomes equal to the 
ice-sublimation temperature, to the radius at which the photo- 
spheric temperature reaches the same value. For example, our 
results indicate that for a standard accretion rate of 1O~ 8 M0 /yr, 
the snow line extends from 0.4 AU (where T e = 160K) to 1.35 
AU (where T m = 160K). We also note that its extent increases 
with mass accretion rate (see also Figure 7). On the other hand, 
for low mass accretion rates (for a 1M Q , 1L Q T Tauri star when 
M< 10- 10 M Q /yr) we find that the snow line is located in the op- 
tically thin region, with entirely different physical and dynam- 
ical characteristics. The consequence of the transition from an 



optically thick to an optically thin snow line will be discussed 
in more detail in a forthcoming paper. 

Thirdly, we note that the predicted overall surface density 
profile is notably shallower than a Minimum Solar Nebula 
model. Only in the marginally opaque region do we reproduce 
the standard MSN radial power law scaling. This will be dis- 
cussed in more detail in §7. 

7. DISCUSSION OF THE RESULTS AND OBSERVATIONAL 
PROSPECTS 

Our new disk model is based on that proposed by CG97, but 
it contains two significant additions. In our analysis, a poly- 
tropic assumption for the disk in optically thick regions enables 
us to study the largely non-isothermal structure of typical inner 
disks, while retaining the analytical nature of the mathematical 
solutions. We have also incorporated the contribution to heat- 
ing due to viscous dissipation which necessarily accompanies 
mass accretion. In addition, unlike CG97 we do not assume a 
MSN surface density profile. We can therefore calculate self- 
consistently all disk quantities as a function of radius assuming 
steady-state accretion or, as we shall consider in a future in- 
vestigation, solve numerically for the temporal evolution of the 
disk as an initial value problem, using the mass accretion rates 
derived in the model 

In this section, we begin by a direct comparison of our model 
predictions with those of CG97 using their chosen fiducial T 
Tauri star characteristics. Then, we critically discuss the model 
assumptions and finally, the model predictions in the light of 
recent observational and theoretical advances. 

7.1. Comparison with the Chiang & Goldreich model 

We compare our model predictions for the temperature pro- 
file in a protostellar disk around a 0.5M Q , 2.5Rq T Tauri star to 
the formula given by CG97 based on their isothermal model. A 
straightforward comparison is impossible since the mass accre- 
tion rate in their model varies between about 5 x 10" 7 to about 
2 x 10~ 8 M Q /yr (assuming the same value of a t in both mod- 
els), as a result of their assumed surface density profile. On the 
other hand, we need to fix M so we choose a fiducial value of 
M = 10~ 8 M Q /yr for the purpose of the comparison. In addition, 
for this simulation only we change the opacity Ky to the value 
of 400cm 2 /g used by CG97. The results are shown in Figure 6. 

We note first that regardless of the inherent difficulty in com- 
paring models with different assumptions, the predictions for 
the photospheric temperatures do, at a first glance, match fairly 
well. This agreement is expected since the only difference in 
the total energy budget in the optically thick, radiatively heated 
regions of the disk arises from slightly different predicted flar- 
ing angles in the two approaches. However, beyond this overall 
match, we note several points of discrepancy. 

Firstly, the photospheric temperature in the optically thick 
inner nebula in our model is in fact systematically lower than 
theirs by 10-20%, except very near the central star. This differ- 
ence solely arises from the lower flaring angle predicted when 
heating is taken into account in the disk mid-plane. 

Secondly, the position of the transition between optically 
thick and optically thin regions of the disk is closer to the cen- 
tral star in their model than in ours. This mis-match is caused 
by our largely different prediction for the surface density pro- 
file of the disk, which implies that the outer regions have higher 
overall mass content (and optical depth) than in a standard MSN 
model for the mass accretion rate considered. 
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Finally, note that their isothermal model predicts a snow 
line located around 1AU, while taking into account the non- 
isothermal nature of the optically thick region implies that the 
snow line extends as far out as 10 AU for this particular T Tauri 
star. Indeed, at 1AU for instance our mid-plane temperature is 
~ 10 3 K, and exceeds the photospheric temperature by a factor 
of five. Neglecting this effect can lead to a gross mis-calculation 
of the position of the snow-line. 

7.2. The validity of the model assumptions 

The three major assumptions of the proposed model are 
the following: viscous dissipation is modeled by an ad-hoc 
a-prescription, the disk stratification is modeled as a polytrope 
in optically thick regions and finally the dust grains are assumed 
to be uniformly mixed with the gas at all heights and radii. 

Turbulent viscous dissipation affects the disk in two ways: it 
promotes the outward transfer of angular momentum enabling 
gas accretion onto the central star, and generates heat preferen- 
tially in the mid-plane of the disk. For modest accretion rates, 
the viscous contribution to the total energy budget of the disk 
at any particular radius is significant only well within the snow- 
line. However, it does does play a role in driving a small tem- 
perature gradient between the mid-plane and the photosphere in 
optically thick regions well outside of the snow-line. By mod- 
eling the turbulent viscous dissipation using an a-prescription 
we implicitly assume that the mechanism driving turbulence is 
robust, and independent of position within the disk. 

While observations ubiquitously suggest that accretion is an 
ongoing process in protostellar disks (Hartmann 2001), and that 
it must be approximately quasi-steady (otherwise, strong ac- 
cumulation of matter near specific locations in the disk would 
be systematically observed), the origin of turbulence still re- 
mains largely unclear. Turbulent mixing cannot be of convec- 
tive nature only, since heating of the surface layers may sup- 
press convection (Watanabe et al. 1990). Mixing by a magneto- 
rotational instability (Balbus & Hawley 1990) is difficult to sus- 
tain in regions where the ionization fraction of the gas is too 
low (the dead-zones). Although surface layers of the disk are 
always potentially subject to the MRI since they are exposed to 
a significant ionizing flux both from cosmic rays and from the 
central star (Gammie 1996), it is not clear that turbulent mo- 
tions driven high above the disk are strong enough to sustain a 
high a-value all the way down to the disk mid plane. Finally, 
turbulent mixing by a sub-critical instability of the background 
Keplerian shear has also been proposed (Richard & Zahn 1999, 
Dubrulle et al. 2005), as well as shearing instabilities in the 
mid-plane layer caused by the sedimentation of the dust parti- 
cles (Goldreich & Ward 1973, Weidenschilling & Cuzzi, 1993 
Garaud & Lin 2004, Gomez & Ostriker 2005), but while the 
former has not yet been formally proven to exist, the typical 
growth rates of the latter may be too long compared with the dy- 
namical timescale to be effective in providing a reliable source 
of angular momentum mixing. 

To complicate matters, it is likely that the intrinsic disk strat- 
ification (both radially and across the mid plane) may result 
in largely non-uniform turbulent transport properties. Thus we 
recognize that the a- prescription is merely an ad-hoc assump- 
tion, necessary to our simplified treatment, and now discuss to 
which extent the results we obtain are model-dependent. 

In our model, the a-prescription is used to calculate the sur- 
face density profile of disks under the assumption of steady- 
state accretion. Consequently, if a varies significantly with ra- 



dius, the radial profiles obtained in the steady-state accretion 
assumption must be modified accordingly, and the radial scal- 
ing laws obtained should be revised entirely. 

The generation of heat by viscous dissipation is on the other 
hand subtly disguised in our model by the assumption of a poly- 
tropic stratification for the gas/dust mixture. Thus the model 
is indirectly dependent on a through the polytropic index as- 
sumed. 

Finally, the last questionable assumption of the model in- 
volves the fully mixed nature of the dust grains with the gas. 
Grains are fully suspended in the gas only when their stopping 
time is much longer than the turnover time of the turbulent ed- 
dies (which no more than one order of magnitude larger than the 
orbital time). For the smallest grains (typically, 0.1 /im) which 
are responsible for most of the absorption of the stellar light 
within the superheated layer, this approximation is largely sat- 
isfactory at all radii. Thus we do not expect the position and 
structure of the superheated layer to be affected by this simpli- 
fied approximation: significant sedimentation of larger grains 
(mm-size) will leave the superheated layer unaffected. On the 
other hand, the disappearance of the smallest grains resulting 
from the overall growth may result in a change in C(r) = z,;/h 
by a factor of order unity. This possibility should be taken into 
account self-consistently in future models. With respect to the 
mid-plane and radial structure, grain growth results in a change 
in k v and (3, which will affect the radial temperature distribu- 
tion (and therefore the spectral slope) in all regions of the disk. 
Sedimentation doesn't affect the re-processing of light in the 
optically thin regions, but may have a detectable effect on the 
position of the photosphere (and therefore the radial profile) in 
the optically thin regions. 

7.3. Discussion of the results and their observational 
implications 

Disk geometry: 

As in CH97, we find that the stellar radiation is primarily ab- 
sorbed by a flaring layer (a > 0) of superheated grains in the 
radiatively heated regions, including in the newly studied case 
of optically thick, non isothermal disks. Regions where viscous 
dissipation is the principal source of energy are, on the other 
hand, found to be non-flaring (a < 0) although their overall ge- 
ometry is very close to flat. 

From the central regions to the outer regions of the disk we 
find that the flaring angle gradually increases outward. The vis- 
cously heated, and radiatively heated, strongly opaque regions 
of the disk are close to a flat geometry (a = 0), with the tempera- 
ture index q increasing from -1 . 125 to -0.85 (corresponding to 
a ~ -0.0625/z/r and a ~ 0.075/z/r respectively). Thus we note 
that a flat disk geometry does not necessarily imply any grain 
sedimentation, and it can also straightforwardly be reproduced 
by a appropriately modeled optically thick disk. 

The marginally opaque transition region (from optically thick 
to optically thin disks) is accompanied by a very sharp, local 
increase in the flaring angle with a = 0.5h/r. This transition 
is therefore accompanied by a sharp local increase in the tem- 
perature, reaching values up to 40% larger than in the optically 
thick region immediately within. 

Self-consistent calculations of the height of the superheated 
layer in the optically thin regions shows that the flaring geome- 
try persists all the way out to several hundred AUs, albeit with 
a slightly smaller flaring angle (a ~ 0.3/z/r). 

Self-shadowing of the outermost regions may occur, as dis- 
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FIG. 6. — Direct comparison between our model predictions (thin lines, line-style coding as in Figure 4) and the CG97 model predictions (thick line) for the 
temperature and surface density in a disk around a typical T Tauri stars with mass 0.5Mq, radius 2.5R@ and effective temperature 4000 K. Note in particular the 
lower effective temperature of our model, the much higher mid-plane temperature, and the higher surface density in the outer regions. The predictions for the 
temperature of the superheated grains are naturally identical. 



cussed earlier, when a exceeds a fraction of the total visual 
optical depth, unless other mechanisms truncating the disk oc- 
cur first. This sharp cut off in the reprocessed radiation could 
appear as an outer edge of the disk. 

Disk SEDs and surface density distribution 
The MIR and FIR photons are emitted by both grains at pho- 
tosphere of the inner opaque regions where T e < 200K and by 
the superheated grains with T+ < 200K at much larger radii. 
Around classical T Tauri stars where M <~ 10~ 7 M Q yr -1 , the 
radii where T e and T+ decrease below 200K are 0.7 and 30 AU 
respectively. Around weak-line T Tauri withM <~ 10~ 9 M Q yr" 1 
the radii where they fall below 2QQK are 0.07 and 30 AU re- 
spectively. Although the area filling factor A s ~ a of the su- 
perheated grains with T± < 200K is well below unity, their total 
emitting area is much larger than the opaque inner region with 
T e < 200K. Consequently, the superheated grains contribute to 
a large fraction of the MIR and FIR radiation and the SED in 
this wavelength range is essentially flat (CG97). 

The magnitude of T e of the opaque regions is generally too 
high for the emerging photons to contribute significantly to the 
sub-mm radiation. These photons are emitted by much-cooler 
and much-larger (mm to cm size) grains near the mid plane of 
marginally- transparent and optically-thin regions in the outer 
disk. At very large disk radii, the superheated grains also con- 
tributes to these long wave radiation. 

A robust feature of our model is the prediction of an overall 
disk surface density much shallower than the standard MSN; 
this prediction is associated with the assumption of steady-state 
accretion. More precisely, we find that the opaque inner disk 
regions satisfy roughly E oc r"° 4 and r"° 7 in the viscously and 
radiatively heated regions respectively. The marginally opaque 
disk regions, however, are radially isothermal which implies a 
surface density profile varying as the MSN profile (Hayashi et 
al. 1985) under the steady-state accretion approximation. This 
is particularly interesting since these intermediate regions can 
have a large radial extent (between 100-200 AU and 4-7 AU 



respectively as M decreases from 10" 7 M Q /yr to 10" 9 M Q /yr). 
Finally, in the optically thin regions we find E oc r~ l 08 . Around 
some T Tauri stars, the outer regions have been resolved by 
mm interferometer observations. Our prediction on the surface 
density distribution in this region is very similar to that inferred 
from the mm maps of resolved protostellar disks around T Tauri 
stars (Mundy et al. 2000, Looney et al. 2003). The disk mass 
determination based on the assumption that these regions are 
optically thin (Beckwith & Sargent 1991) is fully justified. 

Snow line position, structure and history 
As seen in Figure 4, the snow line around 1M Q classical T Tauri 
stars shrinks from a few AU to about 0.6AU as M decreases 
from 10" 7 to 10" 10 MQyr _1 , then re-expands to slightly over 2.2 
AU. This interesting non-monotonous variation of the snow line 
position with mass accretion rate (and therefore with time) is il- 
lustrated in Figure 7. 

For relatively high values of the mass accretion rate, we con- 
firm that the snow line is located in the strongly opaque re- 
gion, and therefore has a significant radial extent. For M ~ 
2 x 10~ 7 Mq y -1 the snow line is located near the present-day 
radius of Jupiter; it may have been critical in promoting the 
emergence of the first gas giant planet at this location in the 
solar system (Ida & Lin 2004). 

During the final stage of disk depletion, the snow line moves 
into the optically thin region. For a 1M Q , 1L T Tauri star this 
transition occurs when M drops below 2 x 10~ 10 M Q yr" 1 . Note 
that the sharpness of the transition observed in Fig 7 may be 
an artefact of the simplistic snow line model used here; we will 
study the physics of the snow line and the consequences of the 
transition from optically thick to optically thick in more detail 
in a subsequent paper. 

Most interestingly, we find that the position of the snow line 
drops down to a radius of about 0.6 AU before expanding out 
again. This location is comparable to the inner boundary of the 
habitable zone around solar-type stars (Kasting et al. 1993). 
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FIG. 7. — Variation of the position of the snow line as a function of mass accretion rate for a IMq, ILq T Tauri star. The solid line represents the location where 
the mid-plane temperature T m = 160K whereas the dash-dotted line represents that where the effective temperature T e = 16QK. Between these two locations, ice can 
retain solid phase in a fraction of the disk near its surface. Since a large fraction of the grains are located near the mid plane, along with the turbulent gas, the ice 
line for the mid-plane temperature is more relevant for determining the phase of most water molecules. The slow decrease in the position of the snow line (for T m ) 
as M initially decreases corresponds to the phase in which the snow line is located in the strongly opaque region. The first plateau in this profile arises when the 
snow line moves into the weakly opaque region. Therefore the observed jump is an artifact of the mathematical asymptotics, and should be interpreted instead as a 
smooth flattening of the profile. The second plateau, however, corresponds to the genuine physical transition of the snow line moving from the optically thick to the 
optically thin region as the disk depletes. 



The discrepancy between our results and those obtained by 
Lecar et al. (2006) for instance arises from our self-consistent 
calculation of the flaring angle in the optically thick regions of 
the disk, combined with the assumption of steady-state accre- 
tion. The consequences of this finding are potentially far reach- 
ing. When the disk mass accretion rate has declined to about 
5 x 1O~ 9 M yr" 1 , the snow line crosses the Earth's orbit and 
steadily continues shrinking; at this point, the local total sur- 
face density of the disk is about 400 g/cm 2 , so that in a ring of 
with a width xr hi n around the Earth's orbit one may find about 
2 x 10~ 2 xM ffl of condensed water in the form of icy particles. 
About half this amount in ices is also found around the orbit of 
Venus when the snow line crosses its orbit. This amount largely 
exceeds the current total water content of Earth and Venus, and 
could provide an alternative explanation for their hydration: the 
possibility of a late stage accretion of ice grains onto large pro- 
toplanetary embryos in situ has the main advantage of provid- 
ing a simple deterministic explanation for the measured D/H 
abundance ratio (Lunine 2005) of the oceans' water. 

Finally, we also note that T Tauri stars younger than 1 Myr 
undergo FU Orionis outbursts (Herbig 1977, Hartmann & 
Keyon 1985). During these outbursts, the central luminosity 
increases above 100L Q (Hartmann 2001) and the snow line can 
near-instantly expand well beyond 100AU, then retreat back to 
its initial position in a short period of time (~ 100 yr). This 
episodic relocation of the snow line will clearly be very disrup- 
tive to the growth of any icy body in most regions of the disk 
within the first 10 5 yr of the disk's evolution. This phenomenon 
may be observable through the coexistence of water vapor and 
ice at large radii in very young protostellar disks. 
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APPENDIX 

This Appendix contains some of the mathematical tools and 
supplements used towards the derivation of many of the equa- 
tions from this paper. 

1. Basic asymptotic integration method: 
Suppose we attempt to integrate a function e"^ (x) within the in- 
terval [a,b], when f(x) has a minimum at x = a. Then 

V /W d*~ j" e -m<x-a)f(a) Ax = e -f(a) 

J a JO 
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-f(a) , , s 
( 1 _ p -ib-a)f'(a)\ 
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(75) 

f(a) V ) f'(a) 

in the limit where (b-a)f'(a) 3> 1. Naturally, the curvature 
terms neglected in the expansion of f(x) near x = a also set an- 
other condition on the applicability of the approximation that 
must be checked for self-consistency. 

When b — > +oo, this approximation can be used for instance 
to recover the well-known approximation to the complementary 
error function as x — > +oo: 



erfc(x) = 



e~ u du ~ 



1 e 



(76) 



This approximation is used in the derivation of equation (23). 
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2. Position of the superheated layer 

The derivation of equation (17) is straightforward using this 
method. Indeed, 

/ Z Kv p(z',r')dl=K V f p m (r')e~^^dr' , (77) 
Jo Jo 
so that we can use the above method with 

a = , b = r 



fir')-- 



2h 2 (r') 



-\n(Zp m (r')) , 



(78) 
(79) 



(note that we thereby consider the possibility of Z having a slow 
variation with r). Here, the function f(r') has a minimum at 
r' = b provided the disk is flaring (i.e. provided h/r is an in- 
creasing function of r). In that case, the formula used is 
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instead, yielding equation (17) a few algebraic manipulations 
later. 

The calculation of the exact position of the superheated layer 
in the limit where the disk is optically thick and the derivation 
of equations (51) and (52) is slightly more involved. To begin 
with, let's rewrite the gas density above the photosphere in this 
region as 

p(r,z) = Pe (r)e~ 21 * , (81) 
where p e (r) is the density at the photosphere, 

(82) 



The equivalent calculation and manipulations that were used 
to obtain equation (17) in the optically thin case can be used 
again to yield this time 

ZpthcKve <- h 2 r dr [ r ) + r dlnr + ^ dr { 2h 2 e 

(83) 

In the limit of a strongly opaque disk, z s 3> h e and z e 3> h e 
(in other words, the atmosphere is much colder than the mid- 
plane), so that we can simplify the above by neglecting the 
second term in the RHS compared to the others. Now we use 
equation (43) to simplify this to 



Zp e h e n v e 2/, » = —q^a , 



£e 



(84) 



where a is defined by equation (18). 

On the other hand, the total surface area of grains above the 
superheated layer is given by 



A s = Zp e K V e 






2 
e 

Zs 



so that 



6 2 h e z s 



(85) 



(86) 



on the assumption (which will be proved below) that z s ~ H. 

To prove that z s ~ H in the limit of a strongly opaque disk 
we construct the quantity 

z 2 -z 2 



H 2 ' 

so that equation (84) can be re-written as 

9 2 a 



Zp m e d K V e 2 ' 5 = 



which can be solved for S as 



2e, 



In 



e n+l ZK V p m H\ 2e 



y8 2 



■In 



9 2 \I(d)a0 2 



(87) 



(88) 



(89) 



which shows that 5 is proportional to e, or in other words that 
z s is close to H when z e is close to H. 

3. Radial structure of the disk: 

When solving for the radial structure of the disk, we are always 
integrating ordinary differential equations of the kind 



d (h 
d7 



= B 



for which the general mathematical solution is 



K+- 



1-fl, 

bTi 1 



-Br 



(90) 



(91) 



However, it is possible to set the integrating constant K to zero 
on physical grounds, since it would otherwise correspond to an 
unphysical constant heating source in the disk. Once this is 
done, we simply get 



h = 



1- 



-B 



(92) 



REFERENCES 



Adachi I., Hayashi C, & Nakazawa K., 1976, Prog. Theor. Phys., 56, 1756 

Adams, F. C, Lada, C. J., & Shu, F. H. 1987, ApJ, 312, 788 

Balbus, S. A, & Hawley, J. E, 1990, ApJ, 376, 214 

Beckwith, S.V.W., & Sargent, A.I. , 1991, ApJ, 382, L31 

Beckwith, S.V.W., & Sargent, A.I. , 1996, Nature, 383, 139 

Bell, K.R., 1999, ApJ, 526, 411 

Calvet, N., D'Alessio, P., Hartmann, L., Wilnet, D., Walsh, A., & Sitko, M., 

2002, ApJ568, 1008 
Chiang, E.I., & Goldreich, P., 1997, ApJ, 490, 368 
Clarke, C. J., Gendrin, A., & Sotomayor, M. 2001, MNRAS, 328, 485 
Ciesla, F.J., & Cuzzi, J.N., 2006, Icarus, submitted, astro-ph/05 11372 
D'Alessio, P., Calvet, N., & Hartmann, L. 2001, ApJ, 553, 321 
Dubrulle, B., Marie, L., Normand, C, Richard, D., Hersant, F, & Zahn, J.P. 

2005, a, 429, 1 
Dullemond, CP, & Dominik, C. 2004, a, 417, 159 

Dullemond, CP, Hollenback, D., & Kamp, I., & D'Alessio, P. 2006, in 

Protostars & planets V, in press. 
Gammie, C.F. 1996, ApJ, 457, 355 
Garaud P., & Lin D.N.C. 2004, ApJ, 608, 1050 



Goldreich P., & Ward W. R., 1973, ApJ, 183, 1051 
Gomez, G.C. & Ostriker, E.C 2005, ApJ, 630, 1093 
Greaves, et al. 2005, ApJ, 619, L187 

Hartmann, L. 2001, Accretion Process in Star Formation, Cambridge 

University Press, Cambridge 
Hartmann, L. & Kenyon, S.J. 1985, ApJ, 299, 462 

Hayashi, C. Nagazawa, K., & Nagagawa, Y. 1985, in Protostars and planets II, 

eds D. Black and M. Mathews, U Arizona Press, 1 100 
Herbig, G. H. 1977, ApJ, 217, 693 
Ida, S., & Lin, D. N. C, 2004, ApJ, 604, 388 

Kasting, H.F., Whitmire, D.P., & Reynolds, R.T. 1993, Icarus, 101, 108 
Lada, C. J., Muench, A. A., Luhman, K. L., Allen, L., Hartmann, L., Megeath, 

T., Myers, P., Fazio, G., Wood, K., Muzerolle, J., Rieke, G, Siegler, N., & 

Young, E. 2006, AJ, 131, 1574 
Laplace, PS. de 1796, Exposition du systeme du monde, Paris 
Lecar, M., Podolak, M., Sasselov, D., & Chiang, E., 2006, ApJ, 640, 1115 
Lin, D.N.C, & Papaloizou, J.C.B., 1980, MNRAS, 191, 37 
Lin, D.N.C, & Papaloizou, J.C.B., 1985, in in Protostars and Planets II, ed D. 

Black & M. Mathews (Tucson: Univ. of Arizona Press), 981 



19 



Looney, L. W., Mundy, L.G. & Welch, W.J., 2003, 592, L255 

Limine, J. , 2005, in Meteorites and tge Early Solar System II eds D. Lauretta, 

L.A. Leshin, & H.Y. McSween Jr., (Tucon: University of Arizona Press) 
Lynden-Bell, D., & Pringle, J.E., 1974, MNRAS, 168, 603 
Marcy, G. W., Cochran, W. D., & Mayor, M. 2000, in Protostars and Planets 

IV, ed V. Mannings, A. P. Boss, & S. S. Russell (Tucson: Univ. of Arizona 

Press), 1285 

Mathis J. S., Rumpl W., & Nordsieck K. H., 1977, ApJ, 217, 425 
McCabe C, Duchene G., & Ghez A.M., 2003, ApJ, 588, 113 
Morrill, G.E., & Voelk, H.J., 1984, ApJ, 287, 371 

Mundy, L. G., Looney, L. W., & Welch, W. J. 2000, in Protostars and Planets 
IV, eds V, Mannings, A.P Boss, & S. S. Russell, (Tucson: Univ. of Arizona 
Press), 355 

Quillen, A.C., 2002, AJ, 124, 400 

Richard, D. & Zahn, J.P 1999, a, 347, 734 



Ruden, S.P, & Lin, D.N.C. 1986, ApJ, 308, 883 
Shakura, N.I., & Sunyaev, R.A., 1973, a24, 337 

Shen, Z.-X., Jones, B., Lin, D. N. C, Liu, X.-W., & Li, S.-L. 2005, ApJ, 635, 
608 

Shuping, R. Y., Bally, J., Morris, M. & Throop, H. 2003, ApJ587, 109 

Supulver K. D., & Lin D. N. C, 2000, Icarus 146, 525 

Suttner, G., Yorke, H. W., & Lin, D. N. C. 1999, ApJ, 524, 857 

Takeuchi, T, & Lin, D. N. C, 2002, ApJ, 581, 1344 

Takeuchi, T, & Lin, D. N. C, 2005, ApJ, 623, 482 

Watanabe, S., Nakagawa, Y, Nakazawa, K., 1990, ApJ, 358, 282. 

Weidenschilling S. J., 1977, MNRAS, 180, 57 

Weidenschilling S. J., & Cuzzi, J. N., 1993, Protostars & Planets III, 1031 
Wilden, B.S., Jones, B.F., Lin, D.N.C., & Soderblom, D.R. 2002, AJ, 124, 2799 
Youdin, A. N, & Chiang, E., 2004, ApJ, 601, 1 109 



